Enhancement of non-local exchange near isolated band-crossings in graphene 
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The physics of non-local exchange interactions in graphene sheets is studied within a ^-orbital tight-binding 
model using a Hartree-Fock approximation and Coulomb interactions modified at short distances by lattice 
effects and at large distances by dielectric screening. We use this study to comment on the strong non-locality 
of exchange effects in systems with isolated band-crossings at energies close to the Fermi energy. We also 
discuss the role of lattice scale details of the effective Coulomb interaction in determining whether or not broken 
symmetry states appear at strong interaction strengths, and in determining the character of those states when they 
do appear. 



I. INTRODUCTION 

Graphene sheets are ideal sp 2 -hybridized pure carbon net- 
works, and have attracted attention in recent years because 
of their appealing combination of theoretical simplicity and 
exceptional physical properties. 1-4 Most electronic properties 
of graphene that have been studied experimentally can be 
successfully described in a non-interacting electron picture. 4 
Electron interaction effects are nevertheless clearly mani- 
fested in perpendicular magnetic fields where they lead to 
quantum Hall ferromagnetism 5 and to the fractional Hall 
effect, 6-8 when twc 9 ' 1(> or more 11 layers are stacked in a way 
which leads to flat bands near the Dirac point, and when 
ribbons with zigzag edges 12 are formed. The strongest in- 
teraction effects so far observed in single-layer graphene in 
the absence of a magnetic field is the logarithmic velocity 
correction 13-16 at momenta near the Dirac point, now apparent 
in photoemission measurements 17 and cyclotron mass mea- 
surements in suspended graphene. 18 Recent Monte Carlo sim- 
ulations of zero-field graphene suggest the more interesting 
possibility of a gap opening 19 at the Dirac point for suffi- 
ciently strong interactions, a property that would drastically 
modify electronic properties. Spontaneous gaps have still not 
been detected in single-layer samples, even when suspended 7 
to reduce disorder and dielectric screening and, although their 
appearance cannot be fully ruled out for cleaner suspended 
samples which might become available in the future, likely do 
not occur. 

In this paper we use a 7C-band lattice-model Hartree-Fock 
calculation to show explicitly that the logarithmic velocity en- 
hancement is related to non-local exchange interactions with 
power-law tails. Our calculations provide a numerical esti- 
mate of the cut-off length which appears in the argument of 
the logarithm in the velocity enhancement expression and can- 
not be obtained from continuum model calculations. We also 
use our calculation to study the role that lattice scale physics 
plays in controlling whether or not gapped states can occur 
in single-layer graphene. We show that the appearance of 
gapped states is sensitive to the long-range of the Coulomb 
interaction. By solving self-consistent ^-orbital Hartree-Fock 
equations, we can assess the possibility of realizing topologi- 
cally non-trivial states like those discussed by Raghu et al.,— 
who study an extended Hubbard model with next neighbor in- 
teractions. 



The paper is organized as follows. We start in section II 
by briefly explaining our implementation of Hartree-Fock the- 
ory for a 7T-orbital lattice model. Here we define our model 
Hamiltonian, comment on how we handle complications due 
to the long-range of the Coulomb interaction, and discuss 
some other technical details of our calculations. In section 
III we carry out a detailed study of the power-law non-local 
exchange interactions and the logarithmic velocity enhance- 
ments they produce. In section IV we present a mean-field 
phase diagram which identifies a variety of distinct broken 
symmetry solutions and captures the dependence of the com- 
petition between them on model parameters. Finally we close 
the paper in section V with a discussion of our findings and 
of the general importance of highly non-local exchange inter- 
actions in semi-metals or semiconductors with isolated band 
crossings, or weakly avoided crossings, close to the Fermi 
level. 



II. 7T-ORBITAL HARTREE-FOCK APPROXIMATION 

The simplest tight-binding model for a carbon lattice retains 
one atomic 2p z orbital on each lattice site and couples them 
with nearest neighbor ppn hopping 21 parameters. We use the 
conventions of Ref. [22], choosing a coordinate system in 
which the honeycomb's Bravais lattice has primitive vectors 

3 1 = a(i,o), ^ 2=a (\'^r)' 

where a = 2.46A is the lattice constant of graphene. The re 
ciprocal lattice vectors are then 
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Because nearest-neighbor hopping connects the honeycomb's 
two triangular sublattices, the 2x2 tight-binding band Hamil- 
tonian is purely off-diagonal: 
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where 70 = — 2.6eV is the hopping parameter and the on-site 
energy has been set to zero. The factor 

/(k) - e 'V/^^l+2 e - ,3 V'/2v^ cos ^^ (4) 
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arises from the phase factors of the Bloch wavefunctions on 
neighboring sites. We have neglected remote neighbor hop- 
ping which gives gives rise to electron-hole asymmetry, i.e. 
to k-dependence of the sum of valence and conduction band 
energies. The convention for Bloch basis state phase factors 
which leads to this form of the Hamiltonian is 

(r|U) = VkA « = -^=E e ' k(R ' +T '^ (r-Ri- n)TJ ff (5) 



where r\ a is the spin part of the wavefucntion, T/ is the position 
of sublattice / in the unit cell, and Nk is the number of unit 
cells in the systemt. The label A = (/, a) combines the lattice 
site label I and the spin label a. 



In this basis the Hartree-Fock Hamiltonian is 



V HF = £ U^' 



kXX' 



L( c k'A' c k'A' 
k' x 
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where (dropping the spin index for simplicity) 



U, 
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(6) 



(q) = (Uk'A'|V|k'AkA') = y*i^Kl('i)%'A(i-i)V(|ri-r 2 |)^(rs)Wa'(^) 
= l£ e i(H(^(»/+v))y (| R , + T/ _ R/ _ T/ ,|) 
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We can simplify the two-body Coulomb integrals in Eqs. 
(IT) and dHJ by combining the momentum-space representa- 
tion for the Coulomb interaction (V 11 (q) = V (q) = 2ne 2 / E r q) 
with the atomic orbital form factor /(q) = / dre~ qr \(j> (r)| . 
We use the explicit form 



f(q) = (l-(r qy)/((l + (r o q) 2 ) 4 ) 



(9) 



obtained by Fourier transforming the radial charge distribu- 
tion of a hydrogenic 2p atomic orbital: 



0(r) = 



1 
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(10) 



The choice a () = a a j ^/30A reproduces the the covalent bond 
radius of carbon a = a/(2y/3j. Calculations in bilayer 

graphene suggest that a larger effective radius ao = 3ao/V30 
is a better choice 10 because it accounts crudely for sp 2 bond- 
ing orbital polarization. The two-body Coulomb integrals are 
then given by 

u h = ~Le iG ' ( ^%(|G|)| 2 n|G|) (11) 

A G 

u x(*0 = ^E e ' G -( T '- T ")|/(|q-G|)| 2 y(|q-G|). (12) 

A G 

where G are the reciprocal lattice vectors and A = N^Aq is the 
system area. 



We will also find it useful to consider an alternate model 
for interactions which assigns a value, V e s(r) to the interac- 
tion strength between electrons which depends only on the 
distance between the lattice sites on which they reside. When 
expressed in terms of V e s(r), 
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(13) 
(14) 



For this real space interaction model we use the simple form 



V eS {d) = l/(e r ^a 2 + d 2 ). 



(15) 



Here a a accounts approximately for the reduction of Coulomb 
interaction strength at short distances due to a orbital polar- 
ization and derealization of the 7£ -charge density on each lat- 
tice sitej22 (In this equation energies are in Hartree {e 2 Iclb) 
units and lengths are in units of the Bohr radius og.) In the 
real space model we choose the on-site interaction parameter 
U separately from the longer range tail; U has been variously 
estimated as having values between U ~ 2eV to U ~ 6eV 24 , 
and up to an effective value of U = 9.3 eVM While Coulomb 
interaction energy at the carbon radius length scale is ~ 20eV, 
and an estimate from the first ionization energy and electron 
affinity of a cabon atom gives U = 9.6eV, 26 the effective on- 
site interaction strength is expected to be greatly reduced in 
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the solid state environment because of screening by polariza- 
tion of bound orbitals on nearby carbon atoms. For larger dis- 
tance interactions we have included a factor 1 je r to account 
for dielectric screening, as in the momentum space version of 
the interaction model. The value chosen for £, can be seen as 
an ad-hoc correction for overestimates of exchange interac- 
tions in Hartree-Fock theory. We study a range of values for 
this interaction parameter model but we believe that a value of 
e, ~ 4 is normally appropriate for graphene placed on a dielec- 
tric substrate. Values chosen for e, and U control not only the 
overall strength of the interaction term 27 but also the relative 
strength of onsite and long range parts of the interaction. We 
will show later how this ratio can play a role in selecting the 
broken symmetry solutions which can appear in these models. 

There are two technical difficulties in these calculations, 
one related to the nature of electron-electron interactions and 
one related to the electronic structure of graphene. The long- 
range of the Coulomb interaction creates some numerical dif- 
ficulties, particularly in evaluating the energies of the charge- 
density-wave states discussed below. We have found the accu- 
rate results can be obtained by choosing a cut-off distance for 
the 1/r tail so that the coupled sites are as nearly as possible 
equally distributed between sublattices. The second challenge 
is related to the band crossing at the Dirac point in graphene, 
at which the wavefunctions which enter the construction of the 
exchange potential have a singular dependence on wavevector. 
Accurate calculations require dense £-point sampling near the 
Dirac point, which increases the computational load rapidly in 
Hartree-Fock calculations because of the non-local exchange 
interactions. In an effort to achieve a satisfactory compro- 
mise between computational load and accuracy we exploit the 
hexagonal symmetry inherent in the problem. This allows us 
to limit our calculations to the irreducible wedge with 1 /12th 
of the Brillouin zone area, even thought the additional phase 
factors in the remainder of the zone need still to be prop- 
erly accounted for when we calculate the exchange poten- 
tial. We use denser adaptive fc-point sampling near the Dirac 
cone while keeping a coarser grid in the remainder of the irre- 
ducible wedge as shown in Fig. Q] In this way it is possible to 
achieve good accuracy while maintaining the numerical load 
at a reasonable level. The coarse k point sampling region was 
typically kept to 16 x 16 density while near the Dirac point we 
have chosen for most of our calculations a sampling density 
corresponding to 512 x 512 points in the full Brillouin zone 
and up to 1024 x 1024 density. 



III. NON-LOCAL EXCHANGE AND LOGARITHMIC 
VELOCITY DIVERGENCE NEAR THE DIRAC POINT 

Graphene's Dirac-like low-energy Hamiltonian 1,28 pro- 
vides an easily studied example of isolated band crossings 
near the Fermi level of a solid. The band crossing at 
the two isolated Fermi points introduce singularities in the 
band Hamiltonian with interesting topological 29 characteris- 
tics, and facilitate the application of field-theoretic pertur- 
bative methods,— As we will discuss later, there are some 
close analogies between interaction physics in graphene and 
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FIG. 1: Choice of the irreducible wedge of the primitive cell (left) 
and the adaptive sampling of the fc-points in the vicinity of the Dirac 
point A* used for most of our calculations. The density of k points 
in the dense region shown in the figure corresponds to a sampling 
density of 512 x 512 points. 



in gapless 30 and narrow gap 31 semiconductors. It has long 
been recognized that interaction effects can become promi- 
nent in gapless, semimetal, and narrow gap systems. For ex- 
ample in a semimetal with a small overlap between valence 
and conduction bands, interactions can induce electron-hole 
pairing and turn the solid into an excitonic insulator i 30,32 In fi- 
nite gap semiconductors Wannier-Mott excitons can form due 
to mutual attraction between a hole and an electron. Non-local 
electron exchange interactions play a relevant role in defining 
the band structure of narrow band semiconductors. 33 In gap- 
less semiconductors exchange induced corrections in the dis- 
persion relation are large near the crossing point and it has 
been argued that virtual generation of excitons can lead to 
a dielectric anomaly.22i2i A general study of materials with 
Fermi points has revealed that for linear band crossings, in- 
teractions always introduce a logarithmically diverging ve- 
locity enhancement, 13 whereas instabilities are expected for 
quadratic crossings. 

The marginal Fermi liquid behavior obtained in 3D, 13 and 
in the graphene 2D case 14 is a consequence of non-local ex- 
change interactions, 1 ^^ as we discuss at length below. To 
demonstrate explicitly how these velocity enhancements ap- 
pear in our calculations we examine the Fock term in Eq. (|6) 
expressed in the sublattice representation: 



Vx(k) = 



v x BA 



(k) 
(k) 



(k) 
(k) 



(16) 



The physics is most clearly explained using the real-space in- 
teraction version of our calculations, although the reciprocal- 
space version is more numerically convenient. The diago- 
nal matrix elements are identical by symmetry and can be 
expressed using the real space sum of effective two body 
Coulomb repulsion in Eq. ( fT~4T >. Using the symmetry prop- 



erty that ( c^, A c k i A ) = ( c k'B c k'B / = 1/2 for every value of k' 
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in neutral graphene, we obtain 

1 N K 

V# A (k) = -^I^ k '- k ) L ^V eff (|Lff|) (17) 

k,i,j 

- -jfe!**.(ii#i)-4 « 8 > 

At half-filling, particle-hole symmetry implies that the 
sublattice-diagonal component of the density-matrix is half of 
the full 7T-band density-matrix, and therefore diagonal in lat- 
tice vector. Only the on-site interaction contributes to V% A (k). 
This contribution to the exchange energy is independent of 
momentum and does not contribute to the quasiparticle veloc- 
ity. For the off-diagonal term, on the other hand, we use the 



relation p AB (k') = (4 g c k / A ) = -/ (k') /2 \f (k') | to obtain 



a) 



b) 




c) 




I D - Odd parity 



2D - Vortex 



3D - Hedgehog 



FIG. 2: Illustration of sublattice psedospin dependence on mo- 
mentum for the Dirac-like Hamiltonians. In ID (a) the pseudospin 
changes direction at the Dirac point, in 2D (b) it has a vortex and in 
3D (c) it has a monopole hedgehog structure. In each case the band 
state sublattice pseudospin changes direction upon crossing through 
the Dirac point. 
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Pab (L/j ) V e ff ( | Lfy I ) . (20) 
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The second form for the right hand side expresses the ex- 
change self-energy explicitly in terms of the sublattice off- 
diagonal element of the real-space density matrix: 



Pab 



— V 



'pAB (k') 



(21) 



In momentum space the Dirac band Hamiltonian's sublattice 
off-diagonal density matrix is singular at the Dirac point be- 
cause the valence band sublattice pseudospin state changes at 
the Dirac point. In a ID model this effect leads to a disconti- 
nuity at the Dirac point, in 2D it leads to momentum space 
vortices, and in 3D to hedgehogs, as illustrated in Fig. [2] 
Because the function /(k) vanishes at the Dirac point, the 
inter-sublattice phase jumps along any line passing through it. 
When this singularity is Fourier transformed to real space it 
leads to a slow power lay decay, as illustrated in Fig. [3]for the 
case of graphene, causing the electron exchange interaction to 
be strongly non-local. 

The behavior of the real space tails can be obtained most 
simply from an analysis of the continuum model. We redefine 
the wave vector k such that it represents the momentum mea- 
sured from the Dirac point K. A general three dimensional 
Hamiltonian with linear dispersion at an isolated band cross- 
ing can be described by the Dirac-Weyl Hamiltonian 



cos ( 



H (k) = hv F ok = Tiv F k ( " m Sm ° e 
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(22) 



The 2D case is obtained by setting 9 = n/2 and ID by setting 
<j> =0,n. For r = in the x direction, we obtain the fol- 
lowing result for the contribution to the density-matrix from a 
valley centered at K: 



Pab (r) 



(2n) 2 J\v 



dk e' kr p AB (k) 



k\<k c 



f_ k dk sgn (k) e\p(ikr) 
Jq c dk kJ\ (kr) 

j£ /2 d 9 Jq c dk sin 2 9 k 2 h (kr sin 9 ) 



ID 

2D 
3D 



(24) 



where d is the dimension of the system and 7j (x) is a Bessel 
function of the first kind. In graphene similar contributions are 
made by the two valleys. The dominant contribution to this 
integral at large r will come from the non-oscillatory kr < 1 
region when J\(x) ~ x/2. Inserting this limit into Eq. ( l24t 
and integrating up to k ~ 1/rwe see that p AB (r) ~ r~ d at large 
r, reminiscent of the dimensional dependence in the decay of 
Friedel oscillations. 36 The off-diagonal density matrix in other 
directions differs only by a phase factor. 

The slow power law decay behavior of the off-diagonal 
density matrix in turn leads to a logarithmic divergence in 
V k V^ B (k) evaluated using Eq. ( fl9t or Eq. ©. We can 
obtain an approximate form for the exchange potentials in Eq. 
d20T i by changing the sum over discrete lattice sites to be a 
continuous integral 



where a = (a x ,a y ,a z ) is the Pauli matrix vector, k = 



kl + k], + kj, tan 9 = ^jk 2 x + kf, jk z and tan = k y /k x . The 
density matrix for the occupied states is then given by 



Vf(k) c -i-/j re -' kr p^(r)y eff (|r|) (25) 



P(k) = i 



(1 -COS0) 



- sin 9e~^ 

(1+COS0) 



(23) 



where £1 is the volume of the unit cell. Using polar coordi- 
nates to represent both k and r we evaluate the radial deriva- 



FIG. 3: (Color online) Absolute value of the sublattice off-diagonal 
density matrix defined in Eq. (12 1 b illustrating the overall power law 
decay r~ 2 . Left Panel: Grey scale map representation where we can 
notice a regular anisotropy of the off diagonal density matrix which 
follows the triangular Bravais lattice structure of graphene. The den- 
sity matrix element falls off with a larger power law along certain dis- 
crete directions. Right Panel: Density matrix at discrete lattice vec- 
tors along the directions rj and V2 indicated in the left panel. For di- 
rection rj we notice a periodic dip in the value of PAfi(r) every three 
lattice constants. The density matrix for these lattice vectors has a 
r~ 3 decay law. The slow r~ 2 power law decay reflects the singular 
dependence of valence band wavefunctions on k at the Dirac points. 
The values of the fitting coefficients are c\ = 0.0037, c% = 0.0015 
and d\ = 0.0019 when distances are measured in nm. 



tive of the exchange potential to obtain 
dVf{k) 1 d 



dk 



2£2 dk 



dre-' kr p AB (r)V eff (\r\ 



= — / dr rcos e 
2D. 



— ikr cos ( 



Pab (r)V eff (|r| 



drr" 



In 



ka 



(26) 



where we integrated the angular variables first, identified the 
lattice constant as the lower limit of the approximate contin- 
uous position integral, and k~ 1 as the upper limit to avoid the 
oscillating regime. Note that the space dimension drops out 
of the final result. Similar conclusions can be reached start- 
ing from Eq. ( [19l > and making a multipolar expansion of the 
Coulomb interaction in momentum space. 

In practical calculations both real-space and reciprocal- 
space Hartree Fock calculations for graphene are able to fol- 
low the velocity enhancement only over a limited range of 
momenta, as illustrated in Fig. [4] The real space formula- 
tion used in the present calculation relies on a truncation of 
the electron interaction range at about six lattice constants, as 
detailed in the appendix. This prescription is able to describe 
a large part of the velocity increase due to non-local interac- 
tions, but saturates more quickly than the momentum space 
calculation which fails at small k values due to the discrete- 
ness of the momentum sums used to construct the exchange 
Hamiltonian. 

In Hartree-Fock continuum model calculations, the 
exchange-enhanced velocity is given by 




0.03 0.03 
k (A- 1 ) 

r < k ► m 




0.001 0.01 0.1 

k (A- 1 ) 

k ► r 



FIG. 4: Upper panel: Tight-binding and Hartree-Fock band struc- 
tures near the Dirac point. The right panels blow up the the small 
rectangular regions shown in the left panels. We observe that the mo- 
mentum space Hartree-Fock calculation (HF2) follows the enhance- 
ment to smaller momenta than the real space truncated interaction 
calculation (HF1). Lower panel: E(k)/k versus k close to the Dirac 
point. The momentum space HF2 calculation used 1024 x 1024 k- 
points in the primitive cell. The velocity enhancement saturates in 
both calculations. The dashed straight line is the small k fit obtained 
using Eq. Ml\ with p c = 30. 



(27) 



where Vf = is the band velocity. The logarith- 

mic enhancement term has the prefactor 0£ ee /4, where a ee = 
e 2 /e,-hVF — (c/e r U/r)a, is the effective fine structure con- 
stant, c is the speed of light, and a is the ordinary vacuum fine 
structure constant. Our full Brillouin zone calculation allows 
us to obtain a numerical value for the dimensionless ultravio- 
let cutoff parameter p c in Eq. f [27b. By fitting the numerical 
results we find that p c . = 30 ± 3. 
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IV. BROKEN SYMMETRY SOLUTIONS PHASE 

DIAGRAM 10 



Recent lattice model Monte Carlo studies of interaction ef- 
fects in graphene carried out by Drut and Lahde 19 predicted 
that they would be strong enough in suspended graphene sam- 
ples to induce a CDW broken symmetry state with different 
electron densities on A and B and a gap emerges in the single- 
particle spectrum. This broken symmetry in graphene is anal- 
ogous to those that supply mass to elementary particles in par- 
ticle physics. It now appears clear that these gaps do not oc- 
cur in experimental samples, possibly because of the role of 
lattice scale physics that is not reliably modeled in these sim- 
ulations. Indeed the size of the gaps must be fixed by ultra- 
violet physics because the two-dimensional Dirac model with 
Coulomb interaction does not define a characteristic energy 
scale. The anticipated broken symmetries do occur in both lat- 
tice and continuum mean-field-theory models of single-layer 
graphene, although the interaction strengths at which they oc- 
cur is likely underestimated by mean-field theory. The calcu- 
lations presented in this section demonstrate that the appear- 
ance or absence of these states is sensitive to lattice model de- 
tail, in particular to the value of the on-site interaction strength 
U and the effective dielectric constant e r . Studies of interac- 
tions based on Hubbard models predict antiferromagnetic in- 
sulating states which appear for U > 2.23 |/o| in Hartree-Fock 
mean- field-theory 3738 and for U > (4.5 ±0.5) \y \ in Quan- 
tum Monte Carlo calculations. 39 A gapped spin-liquid state 
appears for U ~ 3.5 |yo|, 40 before the AF state is reached, in 
the latter case. In graphene, however, any attempt to estimate 
the character of the ground state must account for longer range 
interactions ! 25 ! 41 

For the analysis carried out in this section we have used 
the real space formulation of the effective Coulomb interac- 
tions given in Eqs. (fT3TTT3T > that allows a more direct control 
over the value of the onsite repulsion U and the Coulomb in- 
teraction tail. We used a model with finite truncation of the 
interaction range with a cutoff radius of about six lattice con- 
stants. (Some considerations on optimal cutoff choices are 
explained in the appendix.) Fig. |5]shows the mean-field phase 
diagram produced by these calculations in which both spin- 
density-wave (SDW) and charge-density-wave (CDW) broken 
symmetry states appear. The solid line in the middle of the 
paramagnetic region of this figure follows e r ■ U — 10.2838 
eV. Along this line the Hartree mean field forming a charge 
density state with different densities on A and B sublattices 
vanishes. The ordered states which appear above this line are 
spin-density-wave states, which essentially reflect the physics 
expected for Hubbard models on a square lattice. The or- 
dered states which appear below this line are charge-density- 
wave states. For large U and small e r the charge-density-wave 
boundary is close to the the e, U — 10.2838 eV line, indi- 
cating that its location is determined mainly by this simple 
competition between short-range and long-range interactions. 
When this consideration applies, CDW states cannot occur for 
U > 10.2838 eV since e r cannot take a value smaller than 1 . A 
crude estimate of the onsite repulsion from the carbon atomic 
radius is e 2 ja a ~ 20 eV whereas the value of U that can be ob- 
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1 2 3 4 5 6 7 8 

U (eV) 

FIG. 5: Phase diagram showing where spin-density- wave (SDW) 
and charge-density-wave (CDW) broken symmetry solutions appear 
in our model as a function of the interaction parameters U and £,-. 
Strong short distance repulsion (large U) favors SDW states, whereas 
weak short distance interactions and strong Coulomb interactions 
(small £, ) favors CDW states. Below the the solid line in this figure 
the Hartree mean-field interaction energy is lowered by forming a 
CDW state which has different densities on A and B sublattices. The 
CDW state boundary lies below this line because the band energy fa- 
vors uniform densities. The SDW state is a simple antiferromagnet, 
as expected at large U on bipartite lattices. The arrow in the figure 
shows the critical value U = 2.23 |yo| beyond which SDW solutions 
appear for the pure Hubbard model. The shaded regions in the fig- 
ure indicate the parameter values thought to be most appropriate for 
graphene sheets that are suspended and for those that are supported 
by a dielectric substrate. 



tained from the first ionization potential and electron affinity 
of carbon is U ~ 9.6 eVi 2 ^ The actual value will be further re- 
duced when we account for additional screening effects from 
neighboring and onsite a orbitals, but the physically appro- 
priate value is highly uncertain. In our phase diagram CDW 
solutions, which are favored when the longer range part of the 
interaction is strong but the short-range effective repulsion is 
weak, are restricted to values of e r < 2.2 with small enough 
U. We conclude from this sensitivity that it is not possible to 
reliably predict the occurrence or absence of broken symme- 
try states on the basis of continuum model calculations alone. 
The values of e r and U thought to be appropriate based on con- 
siderations explained elsewher e) 24 ' 25 are consistent with the 
absence of broken symmetry states in single-layer graphene 
samples. 

V. DISCUSSION AND CONCLUSIONS 

In the present work we have presented a detailed analysis of 
mean field Hartree-Fock interaction effects in a lattice model 
of single-layer graphene. We first analyzed the velocity renor- 
malization of the band dispersion near the Dirac point at the 
Hartree-Fock level. These calculations demonstrate explic- 
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itly that the velocity enhancement is produced by non-local 
exchange interactions between different graphene sublattices 
and provide a numerical estimate of a dimensionless ultravi- 
olet parameter which cannot be estimated using Dirac con- 
tinuum model calculations. Similar velocity renormalizations 
occur whenever a linear band crossing occurs at the Fermi 
level producing Fermi points. In dimension d the velocity en- 
hancement is associated with a r~ d power law decay in the 
real space density matrix. Large velocity enhancements will 
also occur for similar reasons whenever band gaps are small, 
or show semimetallic behavior when the character of occu- 
pied states varies rapidly on the scale of the Brillouin-zone, 
although in this case they will always remain finite. This type 
of physics is responsible for strong the non-locality of ex- 
change interaction in gapless or small gap semiconductors 30 
with weak avoided crossing of the bands, in the surface states 
of topological insulators 42 or in metallic armchair carbon 
nanotubes. 43 

The velocity enhancements we explore in graphene are 
partially related to the Fermi surface enhancement incor- 
rectly predicted by Hartree-Fock theory when it is applied to 
metals. 36 In that case the enhancement is always suppressed 
by screening. In graphene, however, the density-of-states van- 
ishes at the Fermi level and screening is less effective ! 49 ' 50 A 
random-phase-approximation theory which includes dynamic 
screening also predicts logarithmic enhancement of the veloc- 
ity, but with a slightly modified logarithm prefactor. 

Our mean field study of broken symmetry states is summa- 
rized by the phase diagram as a function of Coulomb interac- 
tion parameters in Fig. [5] We have shown that CDW states 
are favored by weak on-site interactions and or SDW states 
by strong on-site interactions, but that neither instability oc- 
curs in a broad range of interaction parameter space. The 
most realistic values for the two parameters are still not ac- 
curately known, but may be guessed from the character of 
the broken symmetry states which do in fact occur in the 
quantum Hall regime of graphene in which the kinetic en- 
ergy is quenched 44 . Our suggested values for these parame- 
ters, both for suspended and unsuspended samples are shown 
in Fig. [5] According to the phase diagram we have ob- 
tained, suspended samples of graphene without substrate di- 
electric screening (e r ~ 1) is likely reasonably close to a CDW 
instability. This result is in rough agreement with the lat- 
tice Monte-Carlo calculations of Drut-Lahde 19 who predict a 
band-gap opening for graphene for a critical value of e r ~ 1. 
However, the latest available transport measurements for sus- 
pended graphene- 7 - find a finite resistivity of about 16k£l in 
agreement with early predictions 45 for the minimum conduc- 
tivity for graphene. There is no experimental evidence for an 
insulating CDW state. This discrepancy between experiment 
and present theory signals in part the limitations of 7C-band 
only models that do not include screening of the bare electron 
by carbon a band polarization. An increase of the effective 
dielectric constant from e r = 1 to e r ~ 2 to account for screen- 
ing by degrees-of-freedom not included in the 7T-band model 
would be sufficient to explain the absence of broken symmetry 
states in suspended samples. Recent inelastic X-ray scatter- 
ing experiments 46 in graphite find screening at high energies 



within graphene sheet. These results motivate further efforts 
to estimate high-energy screening in monolayer graphene. 
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Appendix A: Real space truncation of the Coulomb interaction 

We discuss below the optimum choice for the real-space 
interaction cutoff. Even though the definition of effective 
Coulomb integrals in real space has a physically transparent 
meaning, one important drawback is that the long range of 
the Coulomb repulsion makes sums over lattice sites of Eqs. 
([T3i > and (TBI have slow convergence. A simpler method than 
the more accurate Ewald sum 47 consists in introducing a fi- 
nite spherical truncation of the electron interaction range 48 
as an extended Hubbard model where we incorporate farther 
neighbor contributions in the Coulomb term. For many pur- 
poses this method yields correct enough answers because the 
effective reach of the Coulomb interaction shrinks when the 
positive background charge is taken into account. Because 
of the slower decay in real space of the direct Coulomb term 
compared to the exchange potential the inaccuracy in the elec- 
trostatic energy is usually the largest source of error of this 
truncation method specially when there is no charge neutrality 
within the interaction cutoff range in presence of inhomoge- 
neous density distributions. One way to minimize this error is 
to choose the cutoff range such that the electrostatic energy is 
minimized in presence of an symmetric charge imbalance in 
the A and B sublattices of graphene. In order to evaluate the 




FIG. 6: Real space truncation of the interaction range in graphene 
as illustrated with two different cutoff values of L max ~ 2a and 
Lmax ~ 5a. As we change the value of the cutoff radius L max there 
are oscillations in the relative number of carbon lattice sites A and B 
enclosed within the cutoff distance. 



cutoff for the Coulomb interaction term that minimizes the er- 



8 



ror we express the Hartree energy of a CDW state 

where we use the notation Vy = V e ff (dij) for simplicity where 
dij is the distance between the lattice sites i and j. Let us con- 
sider a charge density transfer of Sn from lattice B to lattice 
A such that the densities are = no + Sn and «g = «o — Sn. 
In that case we obtain 



e h = 2 E Vij(n + 5n)(n + Sn) 



ieAjeA 



- E Vy (no - Sn) (n - Sn) 



ieBjeB 



z ieAjeB 

-z E V 'J ("° - ("° + Sn) 

1 ieB.jeA 



(A2) 
(A3) 
(A4) 
(A5) 



The linear terms in Sn above cancel each other and if we 



neglect a constant shift in the origin the electrostatic energy 
difference per lattice is 



SEn, = 



(Sn 



,2 r 



jeA jeB 



(A6) 



where dy is the distance between lattice sites z and jf, t/ = 
V(da) and / is a fixed label belonging to sublattice A. We 
denote the cutoff dependent direct energy corresponding to 
the long ranged part of the Coulomb interaction as 



7 long 

u DI K Ju fnax 



V;, 



(A7) 



jeA,L,„ 



j£B,L,„ 




FIG. 7: Longer ranged contributions to the Hartree energy 
Eq* 8 (L nlax ) as denned in Eq. iAli which shows a strong cutoff 
distance L max dependent oscillation that converges slowly to the lim- 
iting value represented with the horizontal line whose behavior is 
more clearly shown in the inset. We can notice that 
rather close to the asymptotic limit for certain values of L max . A bet- 
ter estimate of the asymptotic limit can be obtained from the behavior 
of Epf°' lg (Lmax) defined in the text. 



which shows an oscillatory dependence on the cutoff dis- 
tance L max > djj as represented in Fig. [7] This behavior 
poses some caveats in extended Hubbard models with only 
one or two neighbor Coulomb interactions when used for ob- 
taining a phase diagram of broken symmetry states involving 
charge density modulations or comparing results between dif- 
ferent models. We can clearly observe that the above men- 
tioned oscillations slowly converge to a constant for very 
large L max . A better estimate for the asymptotic value in 
the limit L max -> ™ can be obtained from E^l ° ng (L max ) = 
Ei=i*£flT S (Lmaxj) /Nmax averaging the values obtained at 
each discrete i th nearest neighbor shell cutoff, where N max 
is the total number of nearest neighbor shells corresponding 
to the cutoff distance L max . We can observe that for certain 
specific values of L max the quantity E^ g (L mux ) is close to 
Ep" 8 (oo). In table I we represent the values of some of these 
select cutoff distances which are the ones that minimize the 
difference in the number of A and B lattices and therefore 
minimizes the deviation from charge neutrality for a CDW 
state within the cutoff range. In our calculations we have used 
a cutoff just above the value L max — 6.429 la listed in the table. 
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